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Abstract. Both in classical and quantum stochastic control theory a major 
role is played by the filtering equation, which recursively updatcs the information 
state of the system under observation. Unfortunately, the theory is plagued by 
infinitc-dimensionality of the information state which severely limits its practical 
applicability, except in a few select cases (e.g. the linear Gaussian case.) One 
solution proposed in classical filtering theory is that of the projection filter. In 
this scheme, the filter is constrained to evolve in a finite-dimcnsional family of 
densities through orthogonal projection on the tangent space with respect to the 
Fisher mètric. Here we apply this approach to the simple but highly nonlinear 
quantum model of optical phase bistability of a stongly coupled two-level atom in 
an optical cavity. We observe near-optimal performance of the quantum projection 
filter, demonstrating the utility of such an approach. 



PACS numbers: 00.00, 20.00, 42.10 

1. Introduction 

Over the past decade it has become increasingly clear that feedback control of quantum 
systems is essentially a problem of stochastic control theory with partial observations 
[1-3]. In this context, the system and observations are generally modeled as a pair 
of Itó (quantum) stochastic differential equations. It is then the goal of the control 
engincer to find a feedback control policy, depending on the system state only through 
the past history of the observations, that achieves a particular control objective. 

In the case of linear system dynamics and observations and Gaussian initial 
conditions, the so-called optimal control problem can be solved exactly both classically 
[4] and quantum-mechanically [1,5] provided that a quadratic performance criterion is 
chosen. This means that the control objective is specified as an optimization problem, 
where a certain cost function (the performance criterion) of the system evolution and 
the control signal is to be minimized. The resulting Linear-Quadratic-Gaussian (LQG) 
control is widely used in many technological applications. An important feature of 
LQG theory is its separation structure: the optimal controller splits up into a part that 
updates the optimal estimate of the system state given the observations (the Kalman 
filter), and an optimal feedback law which is only a function of the state estimate. 

It was originally suggested by Mortensen [6] that the separation structure of LQG 
control carries over even to the nonlinear case. The problem now separates into the 
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nonlinear filtering problcm of finding the optimal estimate of the system statistics 
given the observations and the optimal control problem of finding a feedback law, 
based only on the filtered estimate, that minimizes some performance criterion. The 
estimate propagated by the filter is often referred to as the Information state [7] as 
it contains all information of the system possessed by the observer. Unfortunately, 
nonlinear stochastic control is plagued by two rather severe problcms. First, the 
information state is generally infinite-dimensional even for very simple nonlinear 
systems [8]. Second, even in the finite-dimensional case the nonlinear optimal control 
problem is generally intractable. The latter can sometimes be alleviated by posing a 
less stringent control objective [3]. Nonetheless nonlinear stochastic control remains 
an extremely challcnging tòpic, both in the classical and quantum mechanical case. 

This paper is concerned with the first problcm, that of infmite-dimensionality 
of the nonlinear information state. There is no universal solution to this problcm. 
The most common (though rather ad hoc) approach used by engineers is known as 
the extended Kalman filter [9]. In this scheme, the system dynamics is linearized 
around the current expected system state, and a Kalman filter bascd on the linear 
approximation is used to propagate the estimate. However, aside from the fact that 
the method only performs well for nearly linear systems, it is not clear how it can be 
applied to quantum models§. 

A much more flexible approximation for nonlinear filtering cquations was proposcd 
by Brigo, Hanzon and LeGland [10-12], based on the differential geomètric methods 
of information geometry [13]. In this scheme we fix a finitc-dimcnsional family of 
densities that are assumcd to be good approximations to the information state. Using 
geomètric methods the filter is thcn constrained to evolve in this family. The finite- 
dimensional approximate filter obtained in this way is known as a projection filter, 
and often performs extremely well when the approximating family is chosen wisely. 
Moreover, as this approximate filter is based on the optimal nonlinear filter, instead 
of on the trajectòries of the system state in phase space, it is readily extended to 
the quantum case. Though by no means a universal solution to the filtering problcm, 
we believe that the flexibility and performance of this method likely make it widely 
applicable in the realistic (real-time) implcmcntation of quantum filtering theory. 

In this paper we apply the projection filtering method to a simple, but highly 
nonlinear quantum system: a stongly driven, strongly couplcd two-level atom in a 
resonant single-mode optical cavity [14, 15]. The output field of such an experiment 
exhibits a randomly switching phase, caused by the atòmic spontaneous emission. 
The formalism developed by Brigo et al. can be applied directly to this system if the 
information state (the conditional density of the atom and cavity mode) is represented 
as a Q-function [16]. Remarkably, our projection filter shows strong connections to the 
classical problem of filtering a random jump process in additive white noise [17, 18]. 

Rather than using a quasiprobability representation, a fully quantum theory 
of projection filtering is expressed in terms of finite-dimensional families of density 
operators and quantum information geometry [13]. We will present the general theory 
in a future publication. Nonetheless there is no theoretical objection to the approach 
taken in this paper. In fact, we observo numcrically that the projection filter for our 
model has ncar-optimal performance, demonstrating the utility of this approach. 

§ If the system dynamics can be meaningfully expressed in terms of conjugate pairs of observables, 
one could imagine locally lincarizing the system Langevin equations to obtain a quantum extended 
Kalman filter. To our knowledge this has not yet been attempted. However, it is not clear how to do 
this in c.g. atòmic systems, where the internal degrees of freedom do not obey CCR. 
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Figure 1. Cartoon drawing of the projection filter. (a) The infinite-dimensional space of all densities 
is represented by M, while S is a finite-dimensional submanifold. The filter, defined by a (stochastic) 
difïerential equation in M, flows along a (random) vector field on M. The flow is such that the density 
remains close to S, but S is not invariant. (b) To each point g S the flow associates a (random) 
tangent vector X[0] which has components in both TgS and its complement. The projection filter is 
generated by the vector field in TS that, for each 8, is the orthogonal projection of X[8] onto TgS. 



This paper is organized as follows. In section[2]we introduce the projection filter 
and the neccessary elements of information geometry. Next, in section|3 we introduce 
the physical model that we will be using as an example and obtain the associated 
filtering equation. In section^we obtain the projection filter for our model. Finally, 
in section[3 we present and discuss the results of numerical simulations. 



2. Information geometry and the projection filter 

2.1. The bàsic principle of the projection filter 

The bàsic idea behind the projection filter is illustrated in Figure^ First, we assume 
that the information state can be represented as a probability density, i.e. an integrable 
nonnegative function on some underlying phase space. Though this is not always 
the case even in classical probability, this is generally a good assumption for any 
"reasonable" model. In this paper we will use a well-known quantum quasiprobability 
distribution, the Q-function [16], for this purpose. The set of all possible densities 
forms an infinite-dimensional function space which we will denote by M. 

We also suppose that the information state is well approximated by densities in 
some finite-dimensional subspace S of M. We will assume that S can be given the 
structure of a differential manifold, but we do not require it to be a linear space. As 
such we must be careful in what follows in distinguishing between points in S 1 , points 
in the tangent bundle TS, etc, as in any differential geomètric situation. 

In general, S will not be an invariant set of the filter; if we start with a density in 
S, the filter will cause the density to evolve into a neighborhood of S in M. The idea 
behind the projection filter is simply to constrain the optimal filter to remain in S. 
As S is finite-dimensional, we can then express the projection filter as a differential 
equation in a finite set of local coordinates on S. 

The optimal filter that propagates the information state of the system is given 
by a stochastic differential equation (SDE) in M, and we are seeking to express the 
projection filter as an SDE in S. The precise meaning of this statement is a somewhat 
important point which we will return to at the end of this section; for now, we 
can imagine the filter to be an ordinary differential equation that is driven by the 
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observations, as follows: 

Here pt € M is the information state and Yj is the observation made at time t. X, 
then, is an observation-dependent vector field on M. 

To constrain the filter to evolve in S we must only retain the dynamics of Q that 
is parallel to S; dynamics perpendicular to S will move the density into an vmdesired 
region. Mathcmatically, this idea is simply implcmented if we realize that at each 
point 8 E S, X[8;Y] will have components both in the tangent space TgS and in its 
complement TgS 1 - . We can now constrain the vector field by orthogonally projecting 
X[6; Y] onto TgS for every 8 € S. The resulting approximate filter, in which only the 
dynamics that leaves S invariant is retained, is the projection filter [10-12]. 

Before we can flesh out the details of this scheme we must deal with the fact that 
the filter is not given by a differential equation as in (JTJ , but by an SDE of the form 

dp t = A[ Pt }dt + B[ Pt }dY t (2) 

We would like to think of A + BY t as a "stochastic vector field" so that we can directly 
apply the scheme discussed above. The theory of stochastic differential equations on 
manifolds [19,20] telis us that we can in fact do this, as long as we intèrpret J2J) as a 
Stratonovich SDE 

dp t = A[p t ]dt + B[p t }odY t (3) 

This is not surprising as, for example, Itò's rule is incompatible with the requirement 
that the Lie derivative along a vector field is a derivation [21] (in other words, a 
differential geomètric transformation rule can only contain first derivatives, and the 
only stochastic integral with this property is the Stratonovich integral.) Note that 
usually filtering equations are given in the Itò form; hence we must transform to the 
Stratonovich form before we can derive the projection filter. 



2.2. Information geometry 

In order to perform the key step in the above procedure, the orthogonal projection, we 
need an inner product in the tangent space TgS. A differential manifold is not naturally 
endowed with an inner product structure, however, and hence the projection filter is 
not yet well defined. We need to add to the manifold a Riemannian structure [22]. 
In statistics there is a natural way to do this, and the resulting theory is known as 
information geometry [13]. 

There are different ways of introducing this structure, but perhaps the easiest 
treatment is obtained by considering instead of the densities M the space of square 
roots of densities M 1 / 2 . The fact that any density is integrable guarantees that the 
square root of any density is square integrable; hence M 1 / 2 is a subspace of L 2 , the 
space of square integrable functions, and any vector field on M 1 / 2 takes vàlues in L 2 . 

Similarly, we consider the manifold S 1 ^ 2 , which we will explicitly parametrize as 

S 1/2 = {VpÏ^Ï, 0eecR m } (4) 
That is, S 1 / 2 is a finite-dimensional manifold of square roots of densities, parametrized 
by the local coordinates|| é* <E O. As S 1 ^ 2 C A/ 1 / 2 , for any 9 É 9 the tangent space 

|| By writing this, we are assuming that the entire manifold can be covered by a single coordinate 
chart. Without this assumption the description would be more complicated, as then we couldn't 
dcscribc the projection filter using a simple "extrinsic" SDE in R m . Often we can make our manifold 
obey this property simply by removing a few points; we will see an example of this later. 
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TgS 1 / 2 is the linear subspace of L 2 given by 
T e S 1/2 = Span 



86 1 ' ' d9 r > 



C L 2 (5) 



The reason for working with square roots of densities is that this gives a natural inner 
product in the tangent space, which is simply the Standard L 2 -inner product. In 
particular, we can calculate the associated mètric tensor in the basis of 

/ d^j-ff) d^pj~ff) \ f d^pj^ÏÏ) d^pJ^ÏÏ) 1 

[—a*-' ) = j —m oe^ dx = ï wW (6) 

Up to a factor of 1/4, this is the well-known Fisher information matrix g%j{d). 

We are now in the position to define what we mean by orthogonal projection of 
a vector field on M onto TS. At each 0, the orthogonal projection is 

i=l 3=1 \ / 

where we have used the inverse Fisher information matrix g l ° (9) to account for the 
fact that the basis of (J5J is not orthogonal. This is the main result that is needed to 
obtain projection filters. 

2. 3. Orthogonal projection of a Stratonovich filter 

Let us now discuss how to perform orthogonal projection onto a finite-dimensional 
manifold S for the very general form of a filtering equation. We begin by converting 
the equation to the square root form; this gives 

d^p- t = -^—A[ Pt ]dt+-^—B[ Pt ]odY t (8) 
2^p t 2^/pl 

We now constrain the filter to evolve on S 1 ' 2 through orthogonal projection: 

A\p(.,9 t )] Jx , „ B\p(;6 t )] 



This is just a finite-dimensional SDE for the parameters 9t- To convert the expression 
explicitly into this form, note that by the Stratonovich transformation rule 

dV^A)-Y. d -^W lodel (10) 

i * 

Comparing with (J7J) and JÏÏJ, we find that 



,,(.,<>,) • "/- \ P (.,9,) 

where 

AjM*) = f>^)^-^ (12) 

Equations (|ll·|l . I|12|l and © can be used to directly calculate the projection filter for 
a wide range of models. 
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Figure 2. Schematic of the experimental setup that corresponds to our model. A strongly couplcd 
two-level atom in a resonant, single mode cavity is strongly driven by a resonant driving làser, and 
spontaneously emits in all directions. One of the cavity mirrors is leaky, and the atòmic dynamics is 
observed through homodyne detection of the electromagnètic field in this forward mode. 



3. The physical model and the quantum filter 

3.1. The Jaynes-Cummings model in the strong driving limit 

We consider the following physical system, shown in Figure [3 A two-level atom 
is strongly coupled to the mode of a single-mode cavity. The cavity mode and the 
atòmic frequency are resonant. The atom is strongly driven on resonance by a làser, 
and spontaneously emits in all directions. A forward mode of the electromagnètic ficld 
outside the cavity, initially in the vacuum state, scatters off one of the cavity mirrors. 
By making this mirror slightly leaky, we extract information from the system into 
the external field. Homodyne detection of the forward mode then yields information 
about the atom and cavity. 

The goal of this section is to model this physical system as a pair of Itó quantum 
stochastic diffcrcntial cquations, one of which describes the atom-cavity evolution and 
one describing the homodyne observations. To this end, we begin by writing down 
the full Hamiltonian for the system: 

H = H Q + H d + H 3C + H f + H s (13) 

Here Hq is the free Hamiltonian 

H = tüJoà^a + —^-o- z + / <L·hü)(b\(w)bf((j) + b\(w)b s {w)), (14) 
2 Jo 

Hd, is the drive Hamiltonian 

H d = ih{E/2){e iWoi a - e-^V), (15) 
Hjc is t ne well-known Jaynes-Cummings Hamiltonian 

Hjc = ihg(a'a — acr), (16) 

and Hf and H s are the dipole couplings to the forward and spontaneous emission field 
modes outside the cavity (see e.g. [23]) 

/>oo 

H f = hJ duK f {uj)[ab f f {Lj)+aïb f (Lü)} (17) 

p OC 

H s =h duK s (u)[o-bl(u)+<jïb s {Lü)} (18) 
Jo 

Here a = \g)(e\ is the atòmic lowering operator, a z — [cr,cr\, a is the cavity mode 
lowering operator (we will also use x — a) +a and y — i(a jf — d}), and bf(ui) and b s (w) 
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are the annihilators of the forward and spontaneous emission modes, respectively. The 
resonant frequency of the atom, drive and cavity mode is denoted by uj , £ is the drive 
strength, g is the atom-cavity coupling strength, and Kf(u) and k s (uj) determine the 
frequency-dependent coupling to the external held modes. 

We will assume that lüq £ ^> g > Kf,n s - Following [15], let us switch to the 
interaction picture with respect to Hq + Hd- We obtain the interaction Hamiltonian 

tt rOO 

-i = i s [oV(i) - a<j\t)] + / dw/s.(w)[(j(t)6Í(w)e < (' , '- W0 > t + h.c] 



h 



duj Kf {uj)[ab f f (üj)e l ^-^ n)t + h.c] (19) 



where we have dehncd 

|±)=2- 1 /2(| ff)Tí | e )) ; M =|-)<+|, i*, = \ft,n], (20) 

a(t) = {-i/2)(ne- l£t + l· L z ~ rfé £t ) (21) 

There are two time scales in the Hamiltonian (|19|l . which we will consider separately. 
The first term evolves on the slow time scale of the atòmic evolution. As £ is very 
large compared to the atòmic time scale, we make the rotating wave approximation 
by dropping the rapidly oscillating terms. 

The remaining terms in l|19fl correspond to the fast time scale of interaction with 
the external electromagnètic field. We cannot use the rotating wave approximation for 
these terms, as the external fields are broadband and thus have modes that respond 
on the fast time scale. Instead, we make the weak coupling (Markov) approximation 
for these terms; this results in the following white noise Hamiltonian: 

Hj = H(g/2)fi z x + hV2^ [aÈ\(t) + h.c] + ih[^j+/2 fiÈ f s + (t) 

+ y/^^ÈtJt) + ^/2^Èl_(t) - h.c] (22) 

Here Bp B\ + , È\ z and È\ _ are independent quantum white noises corresponding, 
respectively, to the forward channel and the three spontaneous emission channels at 
lo = lúq + £ , lüq , and ojq — £ (the upper, middle and lower peaks of the Mollow triplet.) 

We refer to [23] for a discussion of the white noise approximation. Care must be 
taken to assign an independent white noise to each frequency component of a(t) (e.g. 
section III. E of [23]); each frequency probes a different subset of the modes b(ui) and 
hence "sees" a different noise. In the weak coupling limit these noises are in fact white 
and independent. For a more rigorous approach to the white noise limit see [24,25]. 
Using the latter approach we can explicitly calculate k — 7tk/(wo) 2 , 7z = ttKs(wo) 2 , 
and 7± = ttk s (u;q ±£ ) 2 . For simplicity, we will assume that approximately 7 Z ,+ - = 7- 

The white noise Hamiltonian (1221) by itself is not well defined. However, we can 
give rigorous meaning to the equation 

dUt i ~ , 

if we intèrpret it as a Stratonovich quantum stochastic diffcrcntial equation [23,25]. 
After conversion to the Itò form, this equation reads 

dU t = [\Fï]2 (/x dBl >+ (t) + n z d,Bl z {t) + fj dBl_(t) -h.c.) 

- iV2~K(a dBj(t) +h.c.) 

- na)adt - (7/2) dt - i(g/2)(j, z xdt] U t (24) 
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Let us now turn to the homodyne observation of the field. The homodyne 
detector measures a quadrature of the forward channel after it has scattered off the 
cavity. We will choose the quadrature B f(t) + B f' (í); the observation process is then 
Y(t) = u}{B f (t) + B f \t))U t (i.e., the photocurrent is I{t) = dY(t)/dt.) Using the 
quantum ítò rules [23,26], we easily find the differential form of this expression: 



We can slightly extend our observation model to account for technical noise, detector 
incfnciency, etc. To model such effects, we add to (|25|) an independent corrupting 
noise cx dC(t) + dC^(t) = dV(t). It is customary in the quantum òptics literature to 
normalizc Y(t) so that dY(t) 2 = dt. In terms of the detection efficiency r\ € (0, 1] 



We will take the Itò equations l|24|) . (|26|l as our model for the system-observation pair. 
3. 2. The quantum filter 

Now that we have a model for the system and the observation process, we can calculatc 
the optimal filter. The derivation of the filtering equation is beyond the scope of this 
paper; for various approaches, see [3,27-30]. We will attempt, however, through 
a simple finite-dimensional analogy, to explain our interpretation of the filtering 
equation, as it is not entirely the same as the interpretation that is often found in 
the physics literature (e.g. [31]). 

The optimal filter propagates the information state, which determines our best 
estimate of every system observable given the observations we have made. In our 
model, every system observable can be represented as a self-adjoint operator X that 
lives on the atom-cavity Hilbert space; as we are working in the Heisenberg picture, 
this observable at time t is given by j t (X) = U\XU t . We must now define what we 
mean by an estimate of an observable. 

The idea behind the concept of estimation is that we have made some observation, 
and given the outeome of this observation we wish to make a guess as to the outeome of 
a different observable that we haven't measured. That is, the estimate of an observable 
X given an observation of Y is some function f(Y) whose outeome represents our 
best guess of X. To find the best estimate we must specify some cost function C[f] to 
optimize; the function that minimizes C is then by definition the optimal estimate. 

The most commonly used estimator is one that minimizes the mean-square error 



The observable f(Y) that minimizes this cost is called the conditional expectation 
£{X\Y) of X given Y. We will use the conditional expectation as our information 
state throughout this paper. However, note that if we had chosen a different cost 
we would obtain a different information state and filter. There is nothing inhcrently 
superior about the choice 1)27(1: in fact, it is sometimes advantageous to choose a 
different estimator with e.g. improved robustness properties [7,32]. 

To understand how the conditional expectation relates to familiar notions from 
quantum theory, we will demonstrate the procedure using a pair of finite-dimensional 
observables [33]. Let X and Y be two n-dimensional observables, n < oo, and let Y 
have m distinct eigenvalues y^. Then Y can be decomposed as 




(25) 



dY(t) = y/2Krju}yU t dt + y/rj[dB f (t) + dBUt)} + y/l-r)dV(t) (26) 



C[f] = ({X-f{Y)f) 



(27) 




(28) 
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where is the projection operator onto the eigenspace corresponding to ï/j. Clearly 
any function of Y is a linear combination of Pj, and vice versa. Hence we identify the 
set of all observables that are functions of Y with the span of {Pi}- The conditional 
expectation £(X\Y) is then the element of this set that minimizes the cost <|27ll . 

To find this element, we use the following trick. The expression (X^Y) defines 
an inner product on the set of n x n complex matrices^f. Using this inner product, we 
orthogonally project X onto the linear space spanned by {Pi}- This gives 



It is a well known fact that the orthogonal projection of some vector v onto a linear 
subspace W with respect to any inner product (a, 6) gives the element w G W that 
minimizes the quantity ((v — w), (v — w)) [34]. In our case, this means that PyX 
minimizes ((X — f*(Y))(X — f(Y))). f(Y) is only an observable, however, if / is real, 
in which case we see that PyX is precisely the conditional expectation £(X\Y). Note 
that the orthogonal projection PyX will always be self-adjoint if X and Y commute. 

Remarkably, whcn X and Y commute, the expression (|29[1 is equivalent to the 
traditional projection postulate. To see this, note that if we observe Y — i/í then PyX 
takes the value (PiX) / (Pi) = (P%X Pi) / (Pi) , which is exactly the expectation of X 
with respect to the initial state projected onto the eigenspace of %ji. The situation is 
somewhat ambiguous for noncommuting X and Y, and we will simply refrain from 
defining the conditional expectation £(X\Y) when [X, Y] ^ 0. 

The quantum filter determines the best estimate of every system observable given 
the observations; i.e., it propagates ^t(X) = £(jt(X)\Y(s < i)), where here £(•[•) is a 
proper infinite-dimensional generalization of (|29|l . A crucial point is that jt{X) and 
Y(s) can in fact be shown to commute for all s < t; this is called the nondemolition 
property by Belavkin [27]. Thus we see that, even though we can evidently intèrpret 
the quantum filter in terms of the projection postulate, we do not need to postulate 
anything beyond the Standard formalism of observables and expectations in quantum 
mechanics. Instead, we see that the filter follows naturally from a statistical inference 
procedurc wherein we find the lcast-squares estimate for every system observable given 
the observations. This point of view is very natural in a control-theoretic context. 

We now give the quantum filter for our model 1241) . (I26|) : we refer to [3,27-30] for 
various approaches to deriving this equation. The result is 



where V[c]X = c) Xc- (é cX + Xe* c) /2. The process dW(t) = dY(t) - y/2~Krj ir t (y) dt 
is known as the innovations process; it describes how "surprised" we are by the 
measurement, as it is the difference between the observation dY(t) and our best 
estimate of what we should observe. It can be shown that, as long as the observation 
process Y(t) has the statistics determined by (|24|) and (|26l) . the innovations process 

% We assumc for simplicity that the expectation map (■) is faithful [33]. If this is not the case, then 
the conditional expectation is not unique. However, all versions of £(-|-) are equivalent in the sense 
that the difference between two versions takes nonzero vàlues with zero probability. 




(29) 



dn t (X) = 7r t ((7/2){X>[/z] + V\ii x ] + V[^]}X) dt 



+ 7T t (2k V[a]X) dt + TT t (i(g/2)[fi z x, X]) dt 
+ ^/^[mtia 1 X - Xa) - -K t (y)n t (X)) x 
(dY(t) - y/2iajirt(i 




(30) 
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dW t is a Wiener process. In some sense this reflects the optimality of the filter, as it 
means that the innovation is unbiased. 

Usually the quantum filter l|30|) is written in its density form. To do this, we 
define a random density operator p t such that + TTt(X) = Tr[Xp t ]. We then find 

dp t = -i(g/2)\p z x,p t ] dt + 2KD[a]p t dt 

+ ( 7 /2){2?M + V[p z ] + V[pï]} Pt dt (31) 

+ \j2nri [ipta' - iap t - Tr[p t y]pt](dY(t) - y / 2K,r/Tr[p t y] dt) 

where T>[c]p = cpé — (c^cp + pc f c)/2. This description is morè economical than the 
raw filter l|3U|) . and appears frequently in the physics literature. We have to be careful, 
however, to intèrpret pt as the information state of an observer with access to Y(t), 
and not as the physical state of the system. This point will be important for the 
interpretation of our results. 

We conclude this section with one more filter, the so-called unnormalized filter, 
which is given by the expression 

dpt = -i{g/2)[p z x, p t ] dt + 2nD\a\pt dt 

+ ( 7 /2){X>M + V[p z ] + T>\jJ]}pt dt (32) 
+ í\J2kt\ \pta) — apt] dY(t) 

The information state pt propagated by this filter is not normalized, Tr[/5t] =/= 1. 
However, it is simply related to the normalized information state by pt = pt/Trlpt]- 
The chief advantage of l-il'li is that it is a linear equation, whereas (|31|l is nonlinear in 
p t . This makes (|32|) somewhat easier to manipulate. 

3.3. The Q -filter 

In [15], it was noticed that density operators of the form 

p = \a)(a\ (8 / dyP a (y)\ty/2)( ïy /2\ (33) 

o=± J 

[\iy/2) are coherent states of the cavity mode] form an invariant set of the filtering 
equation (|31J) . Thus, as long as the initial density is within this set, we can represent 
the filtering equations in terms of the pair of real (Glauber-Sudarshan) functions P ± (y) 
on a line. Substituting (|33() into (|32|l yields the unnormalized P-filter 

dP±{y) = ^-[{±g + Ky)P±{y)]dt 
ày 

+ ~\P?{y) - Ptiv)] dt + ^yPt{y) dY(t) (34) 
For our purposes, it is more convenient to work with unnormalized Q-functions 

Q ± (y) = (± ) iy/2\p\± ) iy/2) = J dy' P ± {y') e~^ 2 ''^ (35) 

+ We mean this in the sense of random variables; that is, Tr[Xpt] is a classical random variable 
with the same statistics as the observable nt(X). We have already implied such a correspondence 
by interpreting Y(t) as a classical stochastic process. In general, we can always express a set of 
observables as classical random variables as long as they commute [33]. 
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as these are always guaranteed to be well-behaved densities [16]. We obtain 

r) 

dy l 



dQtiv) = ~[Qt{y) 



Qt(y)}dt+^-[(±g- 



Ky)Qf{y)]dt 



d 



+ 2 K ^Qf(y)dt 



y- 



dy 



Qf(y)dY(t) (36) 



The simplicity of this expression motivates our choice of this system for demonstrating 
the quantum projection filter. 

Rather than using the Q-function for the projection filter, we could work directly 
with the filter l|31|) in density form and apply methods of quantum information 
geometry [13]. However, note that any mètric on a manifold of densities induces 
a mètric on the corresponding manifold of density operators (e.g. [35]). Thus even the 
Q-function projection filter is a true quantum projection filter, as long as we project 
onto a family of Q-functions that correspond to vàlid quantum states. 



3.4. Observing the spontaneous emission 

Until now we have only observed the forward channel; however, at least in principle, we 
could also observe independently the three spontaneous emission channels B s>z , B St ±. 
We would like to identify a spontaneous emission event with the detection of a photon 
in one of these side channels. As such, in this section we discuss the situation wherein 
direct photodetection is performed in each of the spontaneous emission channels, in 
addition to the homodyne detection of the forward channel. 

The analysis in this case is very similar to the one performed in sections I3.1h 
13.21 The system model is still given by i|24[) . Now, in addition to l|26|l. we need to 
introduce three observation processes N z + _ corresponding to photodetection (with 
perfect efficiency) in the three spontaneous emission channels. The details of this 
setup and the associated filtering equations are well known and we will not repeat 
them here (see e.g. [30]). The full (normalized) filtering equation is given by 

dp t = -i(g/2)\p z x,p t ] dt + 2KV[a]p t dt 

+ ( 1 /2){V[p]+V[p z ]+V[^}}p t dt 

+ \J2kï] [ip t a f - iapt - Tr[p t y]p t ](dY(t) - y / 2Krj Ti[p t y] dt) 
+ Ç[p]p t {dN+{t) - (7/2) TrhíVpt] dt) 
+ G[p z ]pt {dN z (t) - (7/2) dt) 

+ G[p)] Pt (dN_ (t) - (7/2) TrhV Pí ] dt) (37) 

where Q[c]p — cpc* /Tr[cp C t] - p. li can be shown that the statistics of the processes 
-^+,2,-(t) is such that they are counting processes with independent jumps and rates 

(7/2) T^MPí]. (7/2) and ( 7 /2)Tr[MM t Pt], respectively. 

We now have two different filters, equations l|31|l and 1)3 70. for the same physical 
system (|24p. To see how they relate, recali that all the filter is propagating is 
an information state. The information state in (|31|l represents the best estimate 
of an observer who only has access to the homodyne measurement in the forward 
channel. The information state in l|37() . however, represents the best estimate of a 
different observer who has access to both the homodyne observation and to direct 
photodetection of the spontaneous emission channels. Neither information state 
represents the physical state of the system; the latter is given by (|24(l . 
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In practice, the frequency-resolved monitoring of spontaneously emitted photons 
is not (yet) experimentally feasible. Hence we would never use the filter ((37(1 in 
an actual experimental situation. On the other hand, we are able to generate 
photocurrents Y(t), N+ Zl -(t) with the correct statistics in a computer simulation. 
It is then interesting to compare the estimate of an observer who has access to all 
photocurrents to the estimate of a realistic observer who only has access to the forward 
channel. In particular, this gives insight into the question asked in [15], 'in what sense 
should we be able to associate observed phase-switching events [in the forward channel] 
with "actual" atòmic decays?' 

The main reason for introducing (|37l) is that it gives us a convenient way to 
perform computer simulations of the photocurrent Y(t). We wish to generate sample 
paths of Y(t), with the correct statistics, in order to compare the performance of 
the optimal filter 1(31(1 with the projection filter that we will derive shortly. Ideally 
we would directly simulate the system evolution 1)24(1: this problem is essentially 
intractable, however. Fortunately, we have already expressed the statistics of the 
photocurrents Y(t), i\/+ iZ) _(í) completely in terms of the information state. Hence 
we can equivalently simulate the photocurrents by simulating ((37(1 according to these 
rules (dW(t) is a Wiener process, N + (t) has rate (7/2) Tr[p^pp t ], etc.) 

Of course, we can also perform such simulations with (51) . However, the 
advantage of 1(3 7|) is that, if we choose 77 = 1, the pure states are an invariant set 
of this filter. We can thus rewrite the equation as a stochastic Schròdinger equation, 
in which we only have to propagate a vector instead of an operator. This is a much 
more efficient numerical procedure, and is frequently used in quantum òptics [36]. For 
our system, the stochastic Schròdinger equation corresponding to 1(37(1 is given by 

d\i>t) = [{-i{g/2)p z x - in(y) t a - nata - («/4)(y)f) dt 

-V2^(ia+ (y) t /2) dW(t) + {n/{p)p)\' 2 - 1) dN+(t) 

+ (n z - l)dN z (t) + (^/(^t)^ _ i)rfiv_(t)] |^> (38) 

where (c)t = (ipt\c\ïpt), and pt = \ipt)(ipt\- We numerically solve this equation in 
a truncated Fock basis for the cavity mode. The homodyne photocurrent 126(1 is 
calculated from the innovation using dY(t) = \J2nr\ (y) t dt + y/rjdWt + y/1 — rj dV t . 

4. The quantum projection filter 

4-1- The finite-dimensional family 

Before we can obtain a projection filter for ((36(1 . we must fix the finite-dimensional 
family of densities to project onto. Note that each density is actually the pair of Q- 
functions Q ± (y), unlike in section[2]where each density was a single function. However, 
we can easily put the problem into this form by making ± an argument of the function, 
i.e. Q ± (y) = Q(y, ±). The square roots of Q-functions form a perfectly reasonable L 2 
space (1? = L 2 (R) ® L 2 (R)) with the inner product 

/oc 
dyQ{ /2 (y, a)Ql /2 (y, a) (39) 
-°° 

In the following we will use the notations Q^{y) and Q(y, ±) interchangeably. 

Numerical simulations of 1(36(1 show that at any time, both Q + (y) and Q~(y) are 
unimodal, roughly bell-shaped densities with an approximately constant width. This 
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suggests that we can attempt to approximate the information state by unnormalized 
density operators of the form 

p = f+l+X+l <8> |t> + /2)(i/í+/2| + i/ _ |-X-| ® \í/J,-/2)(íij,-/2\ (40) 

This corresponds to the bi-Gaussian family of unnormalized Q-functions 



2^ 



exp 



!/* > 



(41) 



We collect the parameters into a vector 6 — fi ,v ), where 6* € 6 = {/i £ 

R, > 0}. Then the family of square roots of densities 

S 1/2 = {y/q(y,±;9), 8eG} 
is a finite-dimensional manifold in L 2 with the tangent space 

|M|M_,.,} 



T e S 1/2 = Span 
and Fisher mètric 



d^q{y,±;9) dy/q(y,±;0) 



de 1 



Calculating the latter explicitly, we obtain the diagonal matrix 
^)=diag{^, i, ÍL, i. 



(42) 
(43) 

(44) 

(45) 



4-2. The projection filter 

We will perforin projection of the unnormalized filtering equation (|36() . as in [18]. We 
begin by converting the equation into the Stratonovich form: 

dQt(y) = ~[Qf{y) - Qt(y)\ dt + |-[(± 3 + «(i - Hv)Qt{v)] dt 

2 dy 
+ 2/s(l - 2ry)— Q±(y) <ft + «r/(2 - i/ 2 )Q^= (y) dí 



+ V 2k77 



„ 9 
dy 



z (y)odY(t) (46) 



We can now use (|llf) and l|12|) to find dynamical equations for the projection filter 
After tedious but straightforward calculations, we obtain 

\1 
.2 

T7 



dW 



dv t = 

d(4 _ 
dt 

dp£_ = 
dt 



ï^yt — p t)~ Kr l Oh f v t \ dt + \f2Ïvq [j, t v t o dY(t) 



KfJ, t 



2 vt 



(Mi -Mt + ) 



2 Vf 



Conversion to the Itò form changes (|47|l and l|48|) to 



2K< nl· 4vtdY{t) 



-\ v t ~"t )dt+ y/2Íwjfi t v t dY(ï) 



(47) 
(48) 

(49) 
(50) 

(51) 
(52) 
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Finally, we rewrite the equations in terms of the normalized parameters \x and 
vt = v t l( v t + v t)- This g ives 

du+ = - 7 (5+ - 1/2) dt + - P+)0x+ - x 

{dY(t) - [M^ t + + MT(1 - ^} (53) 

^t = -9- + |^(M t - - M + ) (54) 
^L = +g- Kl·l - + 1 (/x+ - /i" ) (55) 

Equations 1|53|) - (|55[1 form the projection filter for our model on the family S 1 / 2 . 

Note that equations (|54|l and (|55|l are singular at = or 1. We can trace 
this back to the fact that we have cheated a little in the definition of our family of 
densities. Whcn v + = (or v~ = 0), the map 9 i— > q(y,±;9) is not invertible, as in 
this case any choice of \x + (or leads to the same density. As we have essentially 
inverted this map to obtain the equations I|53[> - (|55[1 for the parameters, we can hardly 
expect these to be well-defined when this map is not invertible. 

Fortunately the points i/+ = and 1 are never reached if we start the filter with 
< ï> + < 1 . Hence we can make the filter well-defined everywhere simply by removing 
the offending points — and v~ = from S 1 / 2 . The map 6 t— > q(y,±;9) is then 
invertible everywhere (in other words, then the manifold is covered by a single chart.) 
Even if we want to consider starting the filter on = or 1 at t = this is not a 
problem; the filter dynamics will cause v + to evolve off the singular point, so that the 
filter is well defined after an arbitrarily small time step [18]. 



4-3. Connection with the Wonham filter 

There is a remarkable connection between the projection filter obtained in the previous 
section and the theory of jump process filtering. This theory goes back to the beautiful 
clàssic paper by Wonham [17], in which the following problem is solved. 

Denote by x(t) a stationary Mar kovian jump process which switches between two 
states a_ and a + with a rate 7/2; i.e., x(t) is a random telegraph signal [37]. Now 
suppose we do not have access to a complete obeservation of x(t), but only to the 
corrupted observation y(t) defined by 

dy(t) = ^j2~Krjx(t) dt + dw{t) (56) 

where dw(t) is a Wiener process. We can now ask, what is our best guess of the 
probability p+(t) that x{i) = a + , given the observations y(s < í)? The answer is 
given in closed form by (a special case of) the Wonham filter: 

dp+(t) = - 7 b+(í) - 1/2] dt+^2^ P+ (t)[l-p + {t)}(a + - o_) x 

{dy{t) - y/2^[a+p+{t) + a_{\ ~ p+{t))]dt} (57) 

But this is exactly ló-it with /if and /i^ replaced by the constants a + and a_. 

Though intuitively appealing, this is in many ways a remarkable result. There 
appear to be no inherent jumps in either the optimal filter (|31f) . or the system- 
observation pair (|24|l . I|2t)|l from which it was obtained. It is true that we can choose to 
observe a jump process in the spontaneous emission channels, as in l|37(l . but we could 
have equally chosen to perform homodyne or heterodyne detection which do not lead 
to jump process observations. Nonetheless Ij53(l emerges naturally from our model, and 
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Figure 3. Comparison between the optimal and projection filters. A typical observation is shown 
with rj = 1, g = 120, n = 40, 7 = 20, and the integration was performed over 25 000 time steps. The 
top row, figures (a) and (c), were calculated using the optimal filter 1361 : the bottom row, (b) and 
(d), were calculated for the same observation process using the projection filter 1531 — 1531 . 

an expression of the same form can even be obtained directly from (|31|l [15]. Evidently 
there is a deep connection between our system and the theory of jump processes. 

As a classical filter, we can intèrpret the projection filter (|53H - I|55|) as an adaptive 
Wonham filter, where the equations l|54(l - (|55|l continually adapt the parameters a + 
and a_ in the Wonham filter (|53[) . A similar structure was observed in [18], where the 
classical problem of changepoint detection (the detection of a single random jump in 
white noise) was treated using the projection filtering approach. 



5. Numerical results 



In this section we present the results of numerical simulations of the various filters. 
Sample paths of the observation process were generated by numerically solving (|38[1 
with a truncated cavity basis of 25 Fock states and the (appropriately truncated) 
initial state \ipo) = |— ) <g> |0). The thus generated observations were then filtered using 
the optimal filter in Q-function form l|36|) . and using the projection filter (|53fl - l|55[> . 

The optimal filter was implemented using a simple finite-difference scheme [38] on 
a grid of 128 equidistant points in the interval y S [—18, 18], with the appropriately 
truncated initial condition corresponding to l^o)- Finally, the projection filter was 
started with the initial condition ï> + = /i + = fjT = 0, where care was taken not to 
propagate /i + until after the first time step. In all simulations, stochastic integration 
was performed using the stochastic Euler method [39]. 

In Figure 01 a typical filtered sample path is shown. The top row was obtained 
from the optimal filter, while the bottom row was obtained using the projection filter. 
The value inferred for both the conditional probability of finding the atom in the |+) 
state (left column) and the conditional expectation of the y-quadrature (right column) 
are nearly identical for the two filters. Evidently the projection filter is an extremely 
good approximation to the optimal filtering equation. 

Next, in Figure 0] we compare the information state of the optimal and 
projection filters to the information state that is based on the additional observation of 
spontancously emitted photons. The latter filter demonstrates the behavior reported 
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Figure 4. A single run of an experiment is simulated. The dashed line (red) corresponds to the 
optimal estimate of an observer who has access to direct photodetection of the spontaneously emittcd 
photons, as well as homodyne detection of the forward channel. The sòlid line (green) is the optimal 
estimate of a diffcrcnt observer, who only has access to the homodyne photocurrent, for the same 
run of the experiment. The dotted line (blue) is the projection filter estimate based only on the 
homodyne photocurrent. All parameters are the same as in Figuro 131 
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Figure 5. Two more sample paths of the filtered estimate of atòmic state, demonstrating missed 
jumps (top) and false jumps (bottom). All parameters are the same as in Figure 



in [14]. Whenever a photon is observed in one of the side peaks of the Mollow triplet, 
the observer infers that the atom has made a jump. The cstimatcd phase of the cavity 
field then exponentially decays to a steady-state level of (y) = ±g/n. 

If we do not measure the spontaneous emission, our best guess of the atòmic state 
still behaves in a jump-like way However, we see that there is a little delay between 
the time that the observer of spontaneous emission thinks the atom has jumped, and 
the time that the homodyne observer comes to the same conclusion. If we identify 
atòmic decay with the spontaneous emission of a photon, we can now give a fairly 
satisfactory answer to the question posed in [15]: 'In what sense should we be able to 
associate observed phase switching events with "actual" atòmic decays' ? It appears 
that an observed phase switch signals that, had we been making such an observation, 
we would likely have seen a spontaneously emitted photon a little while earlier. 

The detection delay is a rather genèric property of the type of filtcring problems we 
are considering [18,40]. Any time we see a large fiuctuation in the observed process, 
the filter has to decide whether this is a large fiuctuation of the noise, or a large 
fiuctuation of the observed system. As is pointed out by Shiryaev [40] in the context 
of changepoint detection, the filter rides a delicate balance between minimizing the 
delay time and minimizing the probability of "false alarms" . Decreasing the number 
of false alarms (by choosing a different filtering cost) would unavoidably increase the 
delay time, and vice versa. In our system, false alarms are missed jumps and false 
jumps; these do occur, as can be seen in Figure |SJ 

If we wish to generally improve the quality of detection we have no alternative 
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Figure 6. A typical filtered sample path with r) = 1, g = 600, k = 200, 7 = 20. The intcgration 
was performed ovcr 100 000 timc steps. The various line types are as in the previous figures. 



than to increase the signal-to-noise level of the observation. In the case of our system, 
we can do this if we increase g and k while keeping their ratio fixed (the analogy can 
be justified from l(5üjl: the signal x has fixed magnitude x — ±g/n, while the signal-to- 
noise ratio ~ g/yfïï.) A simulation with greatly increased signal-to-noise is shown in 
Figure El In this very strong coupling and damping regime, it appears that not much 
more information can be extracted from observation of the spontaneous emission than 
we could have already inferred from the homodyne photocurrent. 



6. Conclusion 



In this paper, we have suggested that the method of projection filtering can be 
very fruitful whcn applicd to quantum filtering thcory. Using a simple model of a 
strongly coupled two-level atom in a cavity we numerically demonstrated near-optimal 
performance of the projection filter, as is evident from Figures |3HH1 We have also 
shown a connection between this model from cavity QED and the classical Wonham 
filter; the projection filter can be interpreted as an adaptive Wonham filter, applied 
to a quantum model. In future work we will develop a "true" quantum formalism for 
projection filtering, using methods from quantum information geometry. 

The reduction of infinite or high-dimensional filters to a tractable set of equations 
is essential if we wish to perform estimation in real time, for example in a feedback 
control loop. In a control-theoretic context, converting a large, complex system into a 
set of simple equations is known as model reduction. Ideally, however, such a procedure 
should yield some bounds on the error of approximation; in our case we have observed 
numerically that the approximation error is very small, but we have no rigorous bounds 
to back up this statement. In classical control theory of linear systems, the method 
of balanced truncation [41] gives a very general method for model reduction with 
guaranteed error bounds. How to do this effectively for nonlinear systems is still an 
open problem, however, both in classical and in quantum thcory. 
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